function [B,bf] = elliptic_bc(dofs,V, T, TE, ET, d,g,cr,varargin)
%        [B,bf] = elliptic_bc(dofs,d,g,cr,varargin)
% this function is used to treat the boundary condition described by g
% only fit for spline method!

% find the active boundary:
bdr = find(ET(:,2)==0);
nbdr = length(bdr);
% compute the auxilary matrix:
% reserve spaces for boundary matrix
dim = max(max(dofs))-min(min(dofs)) + 1;
Idx1 = zeros(dim,1);
Idx2 = Idx1; row = 1;
bf = zeros(dim,1);
for k=1:nbdr
    eg = bdr(k);
    tri = ET(eg,1); % the father of boundary edge k lies at ET(k,1)!
    degree = d(tri);
    % then i can get the edge's global dof;
    v1_loc = find(TE(tri,:)==eg);
    v2_loc = mod(v1_loc,3)+1;
    v3_loc = mod(v2_loc,3)+1;
    v2 = T(tri,v2_loc);  % global node index for boundary edge
    v3 = T(tri,v3_loc);
    if v1_loc == 2  % in this case, the dofs are sorted inversly 
        v1_loc = -2;
    end
    eg_dof_line = dofs(tri,1) - 1 + cr_indices(0,degree,v1_loc,cr);
    J = ((0:d(tri))/d(tri))';    I = 1-J;
     %comput the  boundary vector and assemble to global vector
    bf(row:row+d(tri)) = vdm11(d(tri))\feval(g,I*V(v2,1)+J*V(v3,1),I*V(v2,2)+J*V(v3,2),varargin{:}); 
    Idx1(row:row+d(tri)) = row:row+d(tri);
    Idx2(row:row+d(tri)) = eg_dof_line;
    row = row + d(tri) + 1;
end
row = row - 1;
bf = bf(1:row);
B = sparse(Idx1(1:row),Idx2(1:row),ones(row,1),row,dim);